function v = composeVFields(v1,v2,method)
% returns v = v1 \circ (id + v2)
if nargin < 3
  method = 'linear';
end
v = zeros(size(v1),class(v1));

switch ndims(v1)-1
 case 3,
  h = eyeHField([size(v,2) size(v,3) size(v,4)]);
  v(1,:,:,:) = v2(1,:,:) + ...
      interp3(squeeze(v1(1,:,:,:)),...
	      v2(2,:,:,:)+h(2,:,:,:),...
	      v2(1,:,:,:)+h(1,:,:,:),...
	      v2(3,:,:,:)+h(3,:,:,:),...
	      method);
  v(2,:,:,:) = v2(2,:,:,:) + ...
      interp3(squeeze(v1(2,:,:,:)),...
	      v2(2,:,:,:)+h(2,:,:,:),...
	      v2(1,:,:,:)+h(1,:,:,:),...
	      v2(3,:,:,:)+h(3,:,:,:),...
	      method);
  v(3,:,:,:) = v2(3,:,:,:) + ...
      interp3(squeeze(v1(3,:,:,:)),...
	      v2(2,:,:,:)+h(2,:,:,:),...
	      v2(1,:,:,:)+h(1,:,:,:),...
	      v2(3,:,:,:)+h(3,:,:,:),...
	      method);
 case 2,
  h = eyeHField([size(v,2) size(v,3)]);
  v(1,:,:) = v2(1,:,:) + interp2(squeeze(v1(1,:,:)),...
				 v2(2,:,:)+h(2,:,:),v2(1,:,:)+h(1,:,:),method);
  v(2,:,:) = v2(2,:,:) + interp2(squeeze(v1(2,:,:)),...
				 v2(2,:,:)+h(2,:,:),v2(1,:,:)+h(1,:,:),method);
 otherwise,
  error('fields must be 2 or 3 dimensional');
end
